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Abstract. We describe a coherent network algorithm for detection and reconstruc- 
tion of gravitational wave bursts. The algorithm works for two and more arbitrarily 
aligned detectors and can be used for both all-sky and triggered burst searches. We de- 
scribe the main components of the algorithm, including the time-frequency analysis in 
wavelet domain, construction of the likelihood time-frequency maps, the identification 
and selection of burst events. 
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1. Introduction 

Coherent network analysis is addressing the problem of detection and reconstruction 
of gravitational waves (GW) with networks of detectors. It has been extensively 
studied in the literature [H El El IU El IE] m application to detection of bursts 
signals, which may be produced by numerous gravitational wave sources in the 
Universe El El HQl HH [121 UHl HH H5] . In coherent methods, a statistic is built 
up as a coherent sum over detector responses. In general, it is expected to be more 
optimal (better sensitivity at the same false alarm rate) than the detection statistics 
of the individual detectors that make up the network. Also coherent methods provide 
estimators for the GW waveforms and the source coordinates in the sky. 

The method we present (called coherent WaveBurst) is significantly different from 
the traditional burst detection methods. Unlike coincident methods [161 El HH] > which 
first identify events in individual detectors by using an excess power statistic and than 
require coincidence between detectors, the coherent WaveBurst method combines all 
data streams into one coherent statistic constructed in the framework of the constrained 
maximum likelihood analysis jl]. Such an approach has significant advantages over the 
coincident methods. First, the sensitivity of the method is not limited by the least 
sensitive detector in the network. In the coherent WaveBurst method the detection is 
based on the maximum likelihood ratio statistic which represents the total signal-to-noise 
ratio of the GW signal detected in the network. Second, other coherent statistics, such as 
the null stream and the network correlation coefficient can be constructed to distinguish 
genuine GW signals from the environmental and instrumental artifacts. Finally, the 
source coordinates of the GW waveforms can be reconstructed. 
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2. Coherent analysis 



The coherent WaveBurst pipeline (cWB) uses a method, for a coherent detection and 
reconstruction of burst signals, based on the use of the likelihood ratio functional jl]. 
For a general case of Gaussian quasi-stationary noise it can be written in the wavelet 
(time-frequency) domain as 

K N 

EE 

k=l i,j=l 



c 



(1) 



where K is the number of detectors in the network, w k [i, j] is the sampled detector data 
(time i and frequency j indices run over some time-frequency area of size N) and ^[i, j] 
are the detector responses. Note, we omit the term 1/2 in the conventional definition of 
the likelihood ratio. The detector noise is characterised by its standard deviation a k [i,j], 
which may vary over the time- frequency plane. The detector responses are written in 
the standard notations 



Zk[hj] = F+kh+{i,j] + F xk h x [i,j] 



(2) 



where F +k (9,(p), F xk (9,<p) are the detector antenna patterns (depend upon source 
coordinates 9 and 0) and h + [i,j], h x [i,j] are the two polarisations of the gravitational 
wave signal in the wave frame. Since the detector responses are invariant with respect 
to the rotation around z-axis in the wave frame, the polarization angle is included in the 
definition of h + and h x . The GW waveforms, h + and h x , are found by variation of C. 
The maximum likelihood ratio statistic is obtained by substitution of the solutions into 
the functional C The waveforms in the time domain are reconstructed from the inverse 
wavelet transformation. Below, for convenience we introduce the data vector w[z, j] and 
the antenna pattern vectors f+[?,j] and f x [^j] 



w[i, j] 



r+(x)^Jj 



w K [i,j] 
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(3) 
(4) 



Further in the text we omit the time-frequency indices and replace the sum ^ =1 with 
J2n > w h ere Qtf is the time-frequency area selected for the analysis. 



The likelihood functional (EqjT]) can be written in the form C = C + + £> 



E 
E 

Qtf l 



W 



W 



U)h + --\f + 



2u2 



f x )h x -^|fx|^ 2 



(5) 
(6) 



where the antenna pattern vectors f + and f x are defined in the Dominant Polarisation 
wave Frame (DPF) [4J. In this frame the antenna pattern vectors are orthogonal to 
each other: (f + • f x ) = 0. The estimators of the GW waveforms are the solutions of the 
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equations 

(w f + ) = |f + |X , (7) 
(wf x ) = \f x \ 2 h x . (8) 

Note, the norms |f + | 2 and |f x | 2 characterise the network sensitivity to the h + and h x 
polarisations respectively. 

2.1. Likelihood regulators 

As first shown in jl], there is a specific class of constraints (often called regulators), 
which arise from the way the network responds to a generic GW signal. A classical 
example is a network of aligned detectors where the detector responses are identical. 
In this case the algorithm can be constrained to search for an unknown function £ rather 
than for two GW polarisations h + and h x , which span a much larger parameter space. 
Note, that in this case |f x | 2 = 0, the Equation [8] is ill-conditioned and the solution for 
the h x waveform can not be found. The regulators are important not only for aligned 
detectors, but also for networks of miss-aligned detectors, for example, the LIGO and 
Virgo network [201 EI] ■ Depending on the source location, the network can be much less 
sensitive to the second GW component (|f x | 2 << |f+| 2 ) and the h x waveform may not 
be reconstructed from the noisy data. 

In the coherent WaveBurst analysis we introduce a regulator by changing the norm 
of the f x vector 

|fx! 2 = |fx| 2 + 5, (9) 

where 5 is a parameter. This is equivalent to adding one more dummy detector to the 
network with the antenna patterns f+,K+i = 0, f X: K+i = Vo and zero detector output 
(xk+i = 0). In this case, the regulator preserves the orthogonality of the vectors f + and 
i' x and the maximum likelihood statistic is written as 

l- = £ f^r^ + =£[<»■ e +> 2 + < w ■ <n . do) 



X I 



where the e + and e' x are unit vectors. Depending on the value of the parameter 5 
different statistics can be generated, for example: 

• 5 = 0- standard likelihood, 

• 5 = oo - hard constraint likelihood. 

2.2. Reconstruction of GW waveforms 



The GW waveforms are given by the solutions of the likelihood functional Eq.5|6 For 
the first GW component the solution is 



4 



When the regulator is introduced it affects the solution for the second GW component. 
In this case we look for such a solution, which gives the second term of the likelihood 
statistic L max , when the solution is substituted into the likelihood functional. Namely, 
we solve the equation 

2(w • f x )h x - |f x fhi - = 0. (12) 

Fx I 

Out of two possible solutions the following one is selected 

^ieWwmS • (13) 



x l 



In case of aligned detectors (|/ x | = 0) this equation results in a trivial solution h x = 0. 



3. Data analysis algorithms 

In this section we describe the algorithms used in the coherent WaveBurst pipeline. 
They include: wavelet transformation, conditioning of input data, construction of time 
delay filters, and generation of coherent triggers. 



3.1. Wavelet transformation 

The discrete wavelet transformations (DWT) are applied to discrete data and produce 
discrete wavelet series w[i,j], where j is the scale index and i is a time index. Applied to 
time series, the DWT maps data from time domain to the wavelet domain. All DWTs 
used in cWB have critical sampling when the output data vector has the same size as 
the input data vector. 

Wavelet series give a time-scale representation of data where each wavelet scale can 
be associated with a certain frequency band of the initial time series. Therefore the 
wavelet time-scale spectra can be displayed as a time-frequency (TF) scallogram, where 
the scale is replaced with the central frequency / of the band. The time series sampling 
rate R and the scale number j determine the time resolution Atj(R) at this scale. The 
DWT preserves the time- frequency volume of the data samples, which is equal to 1/2 for 
the input time series. Therefore the frequency resolution Afj is defined as 1/(2 Atj) and 
determines the data bandwidth at the scale j. For optimal localisation of the GW energy 
on the TF plane, the cWB analysis is performed at several time-frequency resolutions. 

The time-frequency resolution defined above should be distinguished from the 
intrinsic time-frequency resolution of the wavelet transformation, which defines the 
spectral leakage between the wavelet sub-bands and depends on the length of the wavelet 
filter. To reduce spectral leakage we use Meyers wavelets for which long filters can be 
easily constructed [22J. As shown in Figure [TJ it allows us much better localization 
of the burst energy on the time-frequency plane than Symlet60 wavelets used for the 
LIGO S2-S4 analysis [23l 124"] . The disadvantage of the Meyer filters is that for the local 
support they have to be truncated. As the result, the Meyer wavelets are approximately 
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Figure 1. Comparison of spectral leakage from the first (low) frequency band 
0— 1024Hz to the high frequency bands for Haar (black), Symlet60 (red) and Meyerl024 
(green) wavelets after three wavelet decomposition steps. 

orthonormal. From the other side the Meyer filters can be constructed so that the 
Parseval identity holds with better then 0.01% accuracy, which is more than adequate 
for the analysis. 

3.2. Linear prediction error filter 

The linear prediction error (LPE) filters are used to remove "predictable" components 
from an input time series. Usually they are constructed and applied in the time domain. 
In this case the output of the LPE filter is a whitened time series. The LPE filters can 
be also used in the wavelet domain. For construction of the LPE filters we follow the 
approach described in [26J. The symmetric LPE filters can be constructed from the 
backward and forward LPE filters by using classical Levinson algorithm, or the split 
lattice algorithm. 

Since each wavelet layer is a time series, rather then applying the LPE filter to a 
time series x(t), one can perform a wavelet decomposition x(t) — > w(f, t) first, and then 
construct and apply the LPE filter F(f) individually to each wavelet layer. A set of 
filters F(f) removes predictable components (like lines) in the wavelet layers producing 
data w'(t). The filtered time series x'{t) can be reconstructed from w'(t) with the inverse 
wavelet transformation. An example PSD of the filtered segment of S4 data is shown 
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in Figure [2j As one can see, when applied in wavelet domain the LPE filter removes 
spectral lines but preserves the power spectral density of the noise floor. 




O 200 400 600 800 1000 1200 1400 

frequency, Hz 



Figure 2. Power spectra of original (black) and LPE filtered (red) noise of the Hanford 
4k detector. 

3. 3. Time delay filters in the wavelet domain 

The likelihood method requires calculation of the inner products (x n (r n ), x m (r m )), where 
the data streams are shifted in time to take into account the GW signal time delay 
between the detectors n and m. The time delay r n — r m in turn, depends on the source 
coordinates 9 and </>. 

In time domain it is straightforward to account for the time delays. However for 
colored detector noise it is preferable to calculate the maximum likelihood statistics in 
the Fourier or wavelet (time-frequency) domains. In the wavelet domain one needs 
to calculate the inner products (w n (r n ),w m (T m )). The delayed amplitudes can be 
calculated from the original amplitudes (before delay) with the help of the time delay 
filter D kl (r) 

w n>m {i,j,r) = ^2D k i(r,j)w n:m (i + k,j + 1), (14) 

kl 

where k and / are the local TF coordinates with respect to the TF location (i,j). The 
delay filters are constructed individually for each wavelet layer, which is indicated with 
the index j. 

The construction of time delay filters is related to the decomposition of the sampled 
wavelet functions ^>j(t + r) in the basis of the non-shifted wavelet functions ^fj(t). The 
filter construction procedure can be described in the following steps: 

• create a wavelet series with only one coefficient at the TF location (i,j) set to unity, 
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• apply the inverse wavelet transformation reconstructing *&j(t) in time domain, 

• shift *&j(t) by delay time r and perform wavelet decomposition of ^fj(t + t), 

• the resulting wavelet amplitudes at the TF locations (i + k,j + I) give the delay 
filter coefficients Dfc/(r, j) for the wavelet layer j. 

The length of the filter is determined by the requirement on the acceptable energy 
loss when the time delay filter is applied. The fractional energy loss is 

e K = l-Y,D 2 kh (15) 

K 

where the sum is calculated over the K most significant coefficients. The selected 
coefficients are also described by the list of their relative TF locations (k, I) which 
should be stored along with the filter coefficients Dki- Typically K should be greater 
then 20 to obtain the fractional energy loss less than 1%. 



3-4- Generation of coherent triggers 

A starting point of any burst analysis is the identification of burst events (triggers). 
Respectively, this stage of the burst analysis pipeline is called the event trigger generator 
(ETG). Usually, the ETGs based on the excess power statistics of individual detectors 
are used in the analyses [TBI HH US]. Another example of an ETG is the CorrPower 
algorithm [19], which uses cross-correlation between aligned detector pairs to generate 
the triggers. The likelihood statistic used in the coherent WaveBurst utilizes both the 
excess power and the cross-correlation terms. 



3.4-1- Likelihood time-frequency maps In general the likelihood functional is calculated 
as a sum over the data samples selected for the analysis (see Eqjl]). The number of 
terms in the sum depends on the selected TF area in the wavelet domain. When the 
sum consists of only one term, one can write the likelihood functional for a given TF 
location and point in the sky [|J 

C p (i,j,9,(f)) = |w| 2 - |w - f+h+ - f x h x \ 2 . (16) 

Since the entire likelihood approach is applicable to the functional above, one can solve 
the variation problem and find the maximum likelihood statistics L p (9, <fi). They can be 
maximized over the source coordinates 9 and 0, resulting in the statistics 

L m (i,j) = max^{L p (2, j,9,(j))}. (17) 

Calculated as a function of time and frequency, it gives a likelihood time-frequency 
(LTF) map. Figure [3] shows an example of the LTF map for a segment of the S4 data. 

A single data sample in the map is called the LTF pixel. It is characterized by its 
TF location and by the arrays of wavelet amplitudes Wk(i, j,Tk(9,<p)), which are 
used to construct the likelihood statistics L„. 



| For definition of vectors w, f + , and f x see Eq 3 4 
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Figure 3. Example of the likelihood time-frequency map for a magnetic glitch in the 
S4 LlxHlxH2 data. 



3.4-2. Coherent triggers The statistic L m has a meaning of the maximum possible 
energy detected by the network at a given TF location. By selecting the values of L m 
above some threshold, one can identify groups of the LTF pixels (coherent trigger) on 
the time-frequency plane. A coherent trigger is defined for the entire network, rather 
than for the individual detectors. Therefore, further in the text we reserve a name 
"cluster" for a group of pixels selected in a single detector and refer to a group of the 
LTF pixels clS db coherent or network trigger. 

After the coherent triggers are identified, one has to reconstruct the parameters 
of the GW bursts associated with the triggers, including the reconstruction of the 
source coordinates, the two GW polarisations, the individual detector responses and the 
maximum likelihood statistics of the triggers. The likelihood of reconstructed triggers 
is calculated as 

c c (e,<i>) = Y t c p (i,j,e,<f>) (is) 

u 

where the sum is taken over the LTF pixels in the trigger. The maximum likelihood 
statistic L max is obtained by variation of L c over 9 and <fi. Unlike for L p , which is 
calculated for a single LTF pixel, the L max is calculated simultaneously for all LTF 
pixels forming the coherent trigger. 
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4. Selection of coherent triggers 

When the detector noise is Gaussian and stationary, the maximum likelihood L max is 
the only statistic required for detection and selection of the GW events. In this case 
the false alarm and the false dismissal probabilities are controlled by the threshold on 
£max- The real data however, is contaminated with the instrumental and environmental 
glitches and additional selection cuts should be applied to distinguish genuine GW 
signals [HI EH] . Such selection cuts test the consistency of the reconstructed responses in 
the detectors. In the coherent WaveBurst method the consistency test is based on the 
coherent statistics constructed from the elements of the likelihood and the null matrices. 
The likelihood matrix L nm is obtained from the likelihood quadratic form (see 



EqlO) 



nm nm 

where n and m are the detector indexes. The diagonal (off-diagonal) terms of the matrix 
L mn describe the reconstructed normalized incoherent (coherent) energy. The sum of the 
off-diagonal terms is the coherent energy E^ detected by the network. The coherent 
terms can be used to construct the correlation coefficients: 

^nm i j j = ■ \r^) 

V ^ 'nn-L 'mm 

which represent Pearson's correlation coefficients in the case of aligned detectors. We 
use the coefficients r nm to construct the reduced coherent energy 

^coh ^ L nm \r nm | , n ^ ?7i. (21) 

nm 

which provides one of the most efficient selection cuts for rejection of the incoherent 
background events. 

The null matrix represents the normalized energy of the reconstructed noise 

N n m E nm L nm , (22) 

where E nm is the diagonal matrix of the normalized energy in the detectors: E nn = (x^). 
To distinguish genuine GW signals from the instrumental and invironmental glitches we 
introduce the network correlation coefficients 

_ -^coh _ Ccoh ('OQ'v 

net _ Tt HI? f' Cnet _ ~AT I \ ) 

jVull + |-C/coh| ^ v ull + |e C oh| 

where N n \\ is a sum of all elements of the null matrix, which represents the total 
energy in the null stream. Usually for glitches little coherent energy is detected and 
the reconstructed detector responses are inconsistent with the detector outputs which 
results in the large null energy. Therefore the correlation coefficients C nct and c net can 
be used for a signal consistency test which effectively compares the null energy with 
the coherent energy. This is much safer consistency test than the null stream veto [25] 
where the null energy is compared with the estimated noise energy. In any realistic 
data analysis there is always some residual energy left in the null stream. Therefore for 
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strong gravitational waves the energy of the residual signal can be much larger than the 
noise energy that may result in the false rejection of the GW signals. This is not the 
case for the veto cut based on the C net and c aet . 

5. Summary 

In the paper we discussed how the coherent network algorithms are constructed for 
burst searches. We found it convenient to construct coherent burst searches in the 
time-frequency (wavelet) domain, which requires construction of time delay filters. For 
detection of burst signals we combine output of all detectors into one coherent statistic 
- likelihood, which represents the total signal-to-noise ratio of the signal detected in the 
network. To distinguish genuine GW signals from the instrumental and environmental 
glitches we introduced several coherent statistics constructed from the elements of the 
likelihood and null matrices. We do not discuss the performance of the method in 
this paper, however, numerous studies of the method with different sets of LIGO and 
Virgo data have been performed. It was found, in general, that the method has better 
performance than the burst algorithm used for the published analysis of the LIGO 
data [TBI EH EE] • The results of these studies will be presented in subsequent papers. 
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